Vegetation Index Peak Dynamics

This workflow detects the length of time for NDVI and EVI to go from baseline to peak over the course of the year. Each pixel is classified with a value reflecting the length of time in the year for NDVI and EVI to reach peak greenness.

Notes:

Load libraries

library(raster)
## Loading required package: sp
library(rgdal)
## rgdal: version: 1.1-10, (SVN revision 622)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-3
library(rts)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(stringr)
library(foreach)

Create functions

# Raster time series pulled from working directory + data/year/index
make_rts <- function(year = '2012', index = 'ndvi') {
  # Just in case the year is passed as an integer
  year <- as.character(year)
  # **The data directory must be in the current working directory!
  directory <- paste("data", year, index, "processed", sep='/')
  # Get the GeoTIFF files
  files <- list.files(directory, full.names=TRUE, pattern = "*.tif$")
  # Pull all the files into a raster stack.
  r <- stack(files)
  # Extract dates from file names and cast as Date objects
  dates <- as.Date(str_extract(files, '[0-9]{4}-[0-9]{2}-[0-9]{2}'))
  # Reference the raster stack with their respective dates and collect them
  # together into a time series
  index_ts <- rts(r, dates)
  # Give the raster time series result
  return(index_ts)
}

days_to_max <- function(modists, ts_max, timepoints) {
  start_of_year <- timepoints[1]
  # Create an empty raster
  reclass <- raster(crs = ts_max@crs, # Coordinate Ref Sys
                         ext = ts_max@extent, # Extent
                         resolution = res(ts_max), # Resolution
                         vals=NA) # Fill as NAs
  
  # For each raster in TS...
  foreach(i = 1:length(timepoints)) %do% {
    # How many days in?
    days_since <- as.numeric(timepoints[i] - start_of_year)
    # Subtract the max from current raster in time series
    test <- ts_max - modists[[i]]
    # Redefine the 0 values as days since beginning of year
    reclass[test == 0] <- days_since
  }
  return(reclass)
}

annual_stats_raster <- function(modists, FUN) {
  # We have to use the nested index [[1]] to get the layer itself out of the rts object
  return(apply.yearly(modists, FUN)[[1]])
}

Load raster time series

# EVI and NDVI files for the year are moved in their own directory
# e.g. EVI files for 2012 are in data/2012/evi
# Get the time series
ndvi_ts <- make_rts(year = "2012", index = "ndvi")
evi_ts <- make_rts(year = "2012", index = "evi")

Trace/describe NDVI & EVI over the year

# Get rasters of statistics for each pixel over the course of the year
annual_max_ndvi <- annual_stats_raster(ndvi_ts, max)
annual_sd_ndvi <- annual_stats_raster(ndvi_ts, sd)

# Coefficient of variance calculation
# -- not working...?
#coeff_var <- function(x) sd(x)/mean(x)
#coeff_var_ndvi <- annual_stats_raster(ndvi_ts, coeff_var)

annual_max_evi <- annual_stats_raster(evi_ts, max)
annual_sd_evi <- annual_stats_raster(evi_ts, sd)

Reclassify the raster based on how many days it takes to reach max index value

# These are the timepoints of our series
timepoints <- index(ndvi_ts@time)

reclass_ndvi <- days_to_max(ndvi_ts, annual_max_ndvi, timepoints)
reclass_evi <- days_to_max(evi_ts, annual_max_evi, timepoints)

Plot the time series (multiplot), mean curve, SD results, and reclassified rasters for NDVI and EVI

# Dark violet to green
colors <- rev(rainbow(n = 300, s = 1, v = .7, start = .333, end = .8))
# Time series, multiplot
plot(ndvi_ts,
     col=colors)

plot(evi_ts,
     col=colors)

# Mean of all the rasters in each month, 
plot(x = timepoints,
     y = cellStats(ndvi_ts@raster, mean),
     main = "2012 MODIS Raster Tiles\nMean NDVI Value Over Time",
     type='o',
     xlab = "month",
     ylab = "mean NDVI")

plot(x = timepoints,
     y = cellStats(evi_ts@raster, mean),
     main = "2012 MODIS Raster Tiles\nMean EVI Value Over Time",
     type='o',
     xlab = "month",
     ylab = "mean EVI")

plot(annual_sd_ndvi,
     main = "2012 MODIS NDVI\nStandard Deviation for the Year",
     col = colors)

plot(reclass_ndvi,
     main = "2012 MODIS NDVI\nDays to Reach Maximum NDVI",
     col = colors)

plot(annual_sd_evi,
     main = "2012 MODIS EVI\nStandard Deviation for the Year",
     col = colors)

plot(reclass_evi,
     main = "2012 MODIS EVI\nDays to Reach Maximum EVI",
     col = colors)